*************************************
*** Panel Analysis: SHP *************
*************************************

clear
global OUT "SET_YOUR_OUTPATH"

// this is for data which is later read into R
global OUT2 "SET_YOUR_OUTPATH_FOR_DATAFRAMES"	



*create dummy datasets for regsave

*hybrid models
clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_hhinc.dta", replace

clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_hhinc_dec_pp_sqrt.dta", replace

clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_loginc.dta", replace

clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_hhinc_sat.dta", replace

*bivariate
clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_hhinc_0.dta", replace

clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_hhinc_dec_pp_sqrt_0.dta", replace

clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_loginc_0.dta", replace

clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_hhinc_sat_0.dta", replace


*equivalence testing
clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_hhinc_dec_pp_sqrt_e.dta", replace



*lag models
clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_lag_hhinc.dta", replace

clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_lag_hhinc_dec_pp_sqrt.dta", replace

clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_lag_hhinc_sat.dta", replace

clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_lag_loginc.dta", replace

clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_lag_incdiff2.dta", replace

clear
input coef stderr ci_lower ci_upper N _id
. . . . . .
end
save "$OUT2\SHP\shp_coef_lag_inc_subj_obj.dta", replace





*load data
use "SET_YOUR_INPATH\SHP\Data STATA\SHP-Data-W1-W20-STATA\all_waves\shp_long.dta", clear


xtset idpers year


* remove outliers
replace hhinc=. if hhinc>800000 & hhinc<9999999

*use hhinc with standardized hhinc
egen zhhinc = std(hhinc)
replace hhinc = zhhinc

// use income deciles from 1-10
rename hhinc_dec10 hhinc_dec10_
clonevar hhinc_dec10 = hhinc_dec 
 
rename hhinc_dec_pp10 hhinc_dec_pp10_
clonevar hhinc_dec_pp10 = hhinc_dec_pp 
 
rename hhinc_dec_pp_sqrt10 hhinc_dec_pp_sqrt10_
clonevar hhinc_dec_pp_sqrt10 = hhinc_dec_pp_sqrt 


// gen loginc	
gen loginc = ln(hhinc_pp_sqrt)


*create interaction variables
gen hhincXage = hhinc * age18
gen hhinc_ppXage = hhinc_pp * age18
gen hhinc_pp_sqrtXage = hhinc_pp_sqrt * age18

gen hhinc_decXage = hhinc_dec10 * age18
gen hhinc_dec_ppXage = hhinc_dec_pp10 * age18
gen hhinc_dec_pp_sqrtXage = hhinc_dec_pp_sqrt10 * age18

gen hhinc_satXage = hhinc_sat * age18
gen inc_satXage = inc_sat * age18

gen hhinc_decXage1825 = hhinc_dec10 * age1825
gen hhinc_satXage1825 = hhinc_sat * age1825
gen inc_satXage1825 = inc_sat * age1825

gen hhinc_decXage1835 = hhinc_dec10 * age1835
gen hhinc_satXage1835 = hhinc_sat * age1835
gen inc_satXage1835 = inc_sat * age1835
*



* generate income change/jobloss variables

gen incdiff = hhinc_dec-L.hhinc_dec 

recode incdiff (-9/-1=1 "shock")(0/9=0 "no shock")(else=.), gen(incdiff1)
recode incdiff (-9/-2=1 "shock")(-1 0/9=0 "no shock")(else=.), gen(incdiff2)


// jobloss has only 939 cases on 1
gen jobloss = 0
replace jobloss = . if wstat==.
replace jobloss = 1 if wstat==2 & L.wstat==1	

*create engagement scales via factor analysis
reg vote polint polls_times
gen mysample1 = e(sample)
polychoric vote polint polls_times if mysample1==1

global N = r(sum_w)
matrix r = r(R)
factormat r, n($N) factors(1)
predict f1
sum f1
gen engagement = ((f1 - r(min)) / (r(max)-r(min)))*10


* labels
lab var vote "Vote"
lab var polint "Pol. Interest"
lab var polls_times "No. polls part."
lab var engagement "Political Engagement Scale"


**************************************
***** Analysis ***********************
**************************************


estpost summarize vote polint satdem exeff polls_times trust_gov ///
					hhinc hhinc_pp hhinc_pp_sqrt hhinc_dec10 hhinc_dec_pp10 hhinc_dec_pp_sqrt10 hhinc_sat incdiff2 female age18 unemployed nlf edu_med edu_high  migrant  ///
					hhinc_decXage hhinc_satXage hhinc_decXage1835 hhinc_dec_ppXage hhinc_dec_pp_sqrtXage hhinc_satXage1835
esttab using "$OUT\Tables\SHP\summary_stats_shp.rtf", b(2) cells("count mean(fmt(3)) sd(fmt(3))  min max") replace



* regressions with age (continuous)

// for coef export
foreach var of varlist vote polint polls_times  {
	xthybrid `var' hhinc, clusterid(idpers) full
	eststo m001_`var'
	regsave W__hhinc B__hhinc using "$OUT2\SHP\shp_coef_hhinc_0.dta", ci level(95) autoid append
	xthybrid `var' hhinc female age18 unemployed nlf edu_med edu_high migrant, clusterid(idpers) full
	eststo m01_`var'
	regsave W__hhinc B__hhinc using "$OUT2\SHP\shp_coef_hhinc.dta", ci level(95) autoid append
	
	xthybrid `var' loginc, clusterid(idpers) full
	eststo m002_`var'
	regsave W__loginc B__loginc using "$OUT2\SHP\shp_coef_loginc_0.dta", ci level(95) autoid append
	xthybrid `var' loginc female age18 unemployed nlf edu_med edu_high migrant, clusterid(idpers) full
	eststo m02_`var'
	regsave W__loginc B__loginc using "$OUT2\SHP\shp_coef_loginc.dta", ci level(95) autoid append
		
	xthybrid `var' hhinc_dec_pp_sqrt10, clusterid(idpers) full
	eststo m03_`var'
	regsave W__hhinc_dec_pp_sqrt10 B__hhinc_dec_pp_sqrt10 using "$OUT2\SHP\shp_coef_hhinc_dec_pp_sqrt_0.dta", ci level(95) autoid append
	xthybrid `var' hhinc_dec_pp_sqrt10 female age18 unemployed nlf edu_med edu_high migrant, clusterid(idpers) full
	eststo m3_`var'
	regsave W__hhinc_dec_pp_sqrt10 B__hhinc_dec_pp_sqrt10 using "$OUT2\SHP\shp_coef_hhinc_dec_pp_sqrt.dta", ci level(95) autoid append
	
	xthybrid `var' hhinc_sat, clusterid(idpers) full
	eststo m04_`var'
	regsave W__hhinc_sat B__hhinc_sat using "$OUT2\SHP\shp_coef_hhinc_sat_0.dta", ci level(95) autoid append
	xthybrid `var' hhinc_sat female age18 unemployed nlf edu_med edu_high migrant , clusterid(idpers) full
	eststo m4_`var'
	regsave W__hhinc_sat B__hhinc_sat using "$OUT2\SHP\shp_coef_hhinc_sat.dta", ci level(95) autoid append
}

// for equivalence testing
foreach var of varlist vote  polint polls_times {
	xtreg `var' hhinc_dec_pp_sqrt10 age18 unemployed nlf edu_med edu_high, fe
	regsave hhinc_dec_pp_sqrt10 using "$OUT2\SHP\shp_coef_hhinc_dec_pp_sqrt_e.dta", ci level(90) autoid append
}


// interactions
foreach var of varlist vote polint polls_times  {
	xthybrid `var' hhinc_dec_pp_sqrt10 female age18 unemployed nlf edu_med edu_high migrant hhinc_dec_pp_sqrtXage, clusterid(idpers) full
	eststo m7_`var'
}

// test: objective income on subjective income
foreach var of varlist hhinc hhinc_dec_pp_sqrt10 incdiff2 {
	xthybrid hhinc_sat `var' female age18 unemployed nlf edu_med edu_high migrant, clusterid(idpers) full
	eststo inc_so_`var'
}
esttab inc_so_hhinc inc_so_hhinc_dec_pp_sqrt10 inc_so_incdiff2 ///
	using "$OUT\Tables\SHP\hybrid_all_inc_subj_obj.rtf", b(2) se replace
esttab inc_so_hhinc inc_so_hhinc_dec_pp_sqrt10 inc_so_incdiff2 ///
	using "$OUT\Tables\SHP\hybrid_all_inc_subj_obj.tex", b(2) se replace label nonumbers booktabs alignment(D{.}{.}{3}) scalars("N N")
	

*two tables
esttab m3_vote m7_vote m3_polint m7_polint m3_polls_times m7_polls_times ///
	using "$OUT\Tables\SHP\hybrid_all_incdec_pp_sqrt.rtf", b(3) se replace	
esttab m3_vote m7_vote m3_polint m7_polint m3_polls_times m7_polls_times  ///
	using "$OUT\Tables\SHP\hybrid_all_incdec_pp_sqrt.tex", b(3) se replace label nonumbers booktabs alignment(D{.}{.}{3}) scalars("N N")	///
	varlabels ///
				(R__sex "Female" ///
				R__female "Female" ///
				R__migrant "Migrant" ///
				R__edu_med "Edu: medium" ///
				R__edu_high "Edu: high" ///
				W__edu_med "W: Edu: medium" ///
				W__edu_high "W: Edu: high" ///
				B__edu_med "B: Edu: medium" ///
				B__edu_high "B: Edu: high" ///
				W__age18 "W: Age" ///
				B__age18 "B: Age" ///
				W__unemployed "W: Unemployed" ///
				B__unemployed "B: Unemployed" ///
				W__nlf "W: Not in Labor Force" ///
				B__nlf "B: Not in Labor Force" ///
				W__hhinc "W: Income" ///
				B__hhinc "B: Income" ///
				W__hhincXage "W: Inc. * Age" ///
				B__hhincXage "B: Inc. * Age" ///
				W__hhinc_dec_pp_sqrt10 "W: Income (dec.)" ///
				B__hhinc_dec_pp_sqrt10 "B: Income (dec.)" ///
				W__hhinc_dec_pp_sqrtXage "W: Inc. (dec.) * Age" ///
				B__hhinc_dec_pp_sqrtXage "B: Inc. (dec.) * Age" ///
				W__hhinc_sat "W: Sat. Income" ///
				B__hhinc_sat "B: Sat. Income" ///
				W__hhinc_satXage "W: Sat. Inc. * Age" ///
				B__hhinc_satXage "B: Sat. Inc. * Age" ///
				_cons "Constant")		

	

// engagement scale:
xthybrid engagement hhinc_dec_pp_sqrt10 female age18 unemployed nlf edu_med edu_high migrant, clusterid(idpers) full
eststo m_engage	
	
esttab m_engage ///
	using "$OUT\Tables\SHP\hybrid_engage_eqinc.tex", b(3) se replace label nonumbers booktabs alignment(D{.}{.}{3}) scalars("N N")	///
	varlabels ///
				(R__female "Female" ///
				R__migrant "Migrant" ///
				R__edu_med "Edu: medium" ///
				R__edu_high "Edu: high" ///
				W__edu_med "Edu: medium" ///
				W__edu_high "Edu: high" ///
				B__edu_med "Edu: medium" ///
				B__edu_high "Edu: high" ///
				W__age18 "W: Age" ///
				B__age18 "B: Age" ///
				W__unemployed "W: Unemployed" ///
				B__unemployed "B: Unemployed" ///
				W__nlf "W: Not in Labor Force" ///
				B__nlf "B: Not in Labor Force" ///
				W__hhinc_pp_sqrt10 "W: Income (dec.)" ///
				B__hhinc_pp_sqrt10 "B: Income (dec.)" ///
				_cons "Constant")	
	
	
	
	
* vote as logistic
xthybrid vote hhinc_dec_pp_sqrt10 female age18 unemployed nlf edu_med edu_high migrant, clusterid(idpers) family(binomial) link(logit) full
eststo m5_vote_logit
esttab m5_vote_logit using "$OUT\Tables\SHP\hybrid_vote_logit.rtf", b(3) se replace
esttab m5_vote_logit using "$OUT\Tables\SHP\hybrid_vote_logit.tex", b(3) se replace label nonumbers booktabs alignment(D{.}{.}{3}) scalars("N N")	///
	varlabels ///
				(R__sex "Female" ///
				R__female "Female" ///
				R__migrant "Migrant" ///
				R__edu_med "Edu: medium" ///
				R__edu_high "Edu: high" ///
				W__edu_med "W: Edu: medium" ///
				W__edu_high "W: Edu: high" ///
				B__edu_med "B: Edu: medium" ///
				B__edu_high "B: Edu: high" ///
				W__age18 "W: Age" ///
				B__age18 "B: Age" ///
				W__unemployed "W: Unemployed" ///
				B__unemployed "B: Unemployed" ///
				W__nlf "W: Not in Labor Force" ///
				B__nlf "B: Not in Labor Force" ///
				W__hhinc "W: Income" ///
				B__hhinc "B: Income" ///
				W__hhincXage "W: Inc. * Age" ///
				B__hhincXage "B: Inc. * Age" ///
				W__hhinc_dec_pp_sqrt10 "W: Income (dec.)" ///
				B__hhinc_dec_pp_sqrt10 "B: Income (dec.)" ///
				W__hhinc_dec_pp_sqrtXage "W: Inc. (dec.) * Age" ///
				B__hhinc_dec_pp_sqrtXage "B: Inc. (dec.) * Age" ///
				W__hhinc_sat "W: Sat. Income" ///
				B__hhinc_sat "B: Sat. Income" ///
				W__hhinc_satXage "W: Sat. Inc. * Age" ///
				B__hhinc_satXage "B: Sat. Inc. * Age" ///
				_cons "Constant")








*************************
* FE models
*************************

***** Lagged IVs and Dvs

* restriction: variables need to be at least in 4 waves

	
* FE model
foreach var of varlist vote polint polls_times  {
	xtreg `var' L(0/3).hhinc F1.hhinc age18 unemployed nlf edu_med edu_high, fe
	eststo fe0_`var'_inc
	regsave F1.hhinc hhinc L1.hhinc L2.hhinc L3.hhinc using "$OUT2\SHP\shp_coef_lag_hhinc.dta", ci level(95) autoid append
	xtreg `var' L(0/3).loginc F1.loginc age18 unemployed nlf edu_med edu_high, fe
	eststo fe0_`var'_loginc
	regsave F1.loginc loginc L1.loginc L2.loginc L3.loginc using "$OUT2\SHP\shp_coef_lag_loginc.dta", ci level(95) autoid append
	xtreg `var' L(0/3).hhinc_dec_pp_sqrt10 F1.hhinc_dec_pp_sqrt10 age18 unemployed nlf edu_med edu_high, fe
	eststo fe_`var'_inc_pp_s
	regsave F1.hhinc_dec_pp_sqrt10 hhinc_dec_pp_sqrt10 L1.hhinc_dec_pp_sqrt10 L2.hhinc_dec_pp_sqrt10 L3.hhinc_dec_pp_sqrt10 using "$OUT2\SHP\shp_coef_lag_hhinc_dec_pp_sqrt.dta", ci level(95) autoid append
	xtreg `var' L(0/3).hhinc_sat F1.hhinc_sat age18 unemployed nlf edu_med edu_high, fe
	eststo fe_`var'_inc_sat
	regsave F1.hhinc_sat hhinc_sat L1.hhinc_sat L2.hhinc_sat L3.hhinc_sat using "$OUT2\SHP\shp_coef_lag_hhinc_sat.dta", ci level(95) autoid append
	xtreg `var' L(0/3).incdiff2 F1.incdiff2 age18 unemployed nlf edu_med edu_high, fe
	eststo fe_`var'_incdiff2
	regsave F1.incdiff2 incdiff2 L1.incdiff2 L2.incdiff2 L3.incdiff2 using "$OUT2\SHP\shp_coef_lag_incdiff2.dta", ci level(95) autoid append
}


esttab fe_vote_incdiff2 fe_polint_incdiff2 fe_polls_times_incdiff2  ///
	using "$OUT\Tables\SHP\fe_lag_all_incdiff2.rtf", b(3) se replace	
esttab fe_vote_incdiff2 fe_polint_incdiff2 fe_polls_times_incdiff2 ///
	using "$OUT\Tables\SHP\fe_lag_all_incdiff2.tex", b(3) scalars("N N") se replace label nonumbers booktabs alignment(D{.}{.}{3})	///
	varlabels ///
				(edu_med "Edu: medium" ///
				edu_high "Edu: high" ///
				age18 "Age" ///
				unemployed "Unemployed" ///
				nlf "Not in Labor Force" ///
				hhinc "Income\textsubscript{t0}" ///
				L.hhinc "Income\textsubscript{t-1}" ///
				L2.hhinc "Income\textsubscript{t-2}" ///
				L3.hhinc "Income\textsubscript{t-3}" ///
				F.hhinc "Income\textsubscript{t+1}" ///
				hhinc_dec_pp_sqrt10 "Income (dec.)\textsubscript{t0}" ///
				L.hhinc_dec_pp_sqrt10 "Income (dec.)\textsubscript{t-1}" ///
				L2.hhinc_dec_pp_sqrt10 "Income (dec.)\textsubscript{t-2}" ///
				L3.hhinc_dec_pp_sqrt10 "Income (dec.)\textsubscript{t-3}" ///
				F.hhinc_dec_pp_sqrt10 "Income (dec.)\textsubscript{t+1}" ///
				hhinc_sat "Sat. Income\textsubscript{t0}" ///
				L.hhinc_sat "Sat. Income\textsubscript{t-1}" ///
				L2.hhinc_sat "Sat. Income\textsubscript{t-2}" ///
				L3.hhinc_sat "Sat. Income\textsubscript{t-3}" ///
				F.hhinc_sat "Sat. Income\textsubscript{t+1}" ///
				incdiff2 "Inc. Shock\textsubscript{t0}" ///
				L.incdiff2 "Inc. Shock\textsubscript{t-1}" ///
				L2.incdiff2 "Inc. Shock\textsubscript{t-2}" ///
				L3.incdiff2 "Inc. Shock\textsubscript{t-3}" ///
				F.incdiff2 "Inc. Shock\textsubscript{t+1}" ///
				_cons "Constant")
	
		
*Cross-check: objective income on subjective income 
foreach var of varlist hhinc hhinc_dec_pp_sqrt10 incdiff2 {
	xtreg hhinc_sat L(0/3).`var' F1.`var' age18 unemployed nlf edu_med edu_high, fe
	eststo m_`var'
	*regsave F1.`var' `var' L1.`var' L2.`var' L3.`var' using "$OUT2\SHP\shp_coef_lag_inc_subj_obj_`var'.dta", ci level(95) autoid append
}
esttab m_hhinc m_hhinc_dec_pp_sqrt10 m_incdiff2   ///
	using "$OUT\Tables\SHP\fe_lag_all_inc_subj_obj.rtf", b(3) se replace	
esttab m_hhinc m_hhinc_dec_pp_sqrt10 m_incdiff2   ///
	using "$OUT\Tables\SHP\fe_lag_all_inc_subj_obj.tex", b(3) se replace label nonumbers booktabs alignment(D{.}{.}{3}) scalars("N N")

	
	
* Error correction model: lagged DV	

foreach var of varlist vote polint polls_times  {
	
	// exogenous
	xtdpdsys `var' hhinc_dec_pp_sqrt10 age18 unemployed nlf edu_med edu_high, twostep vce(robust) lags(2) maxldep(2) maxlags(2) artests(2)
	estat abond
	eststo ec_exo_`var'_inc_pp_s
	
	// predetermined
	xtdpdsys `var' age18 unemployed nlf edu_med edu_high, twostep vce(robust) lags(2) maxldep(2) maxlags(2) artests(2) pre(hhinc_dec_pp_sqrt10)
	estat abond
	eststo ec_pre_`var'_inc_pp_s
	
	// endogeneous
	xtdpdsys `var' age18 unemployed nlf edu_med edu_high, twostep vce(robust) lags(2) maxldep(2) maxlags(2) artests(2) end(hhinc_dec_pp_sqrt10)
	estat abond
	eststo ec_end_`var'_inc_pp_s
}


esttab ec_exo_vote_inc_pp_s ec_pre_vote_inc_pp_s ec_end_vote_inc_pp_s ///
	ec_exo_polint_inc_pp_s ec_pre_polint_inc_pp_s ec_end_polint_inc_pp_s ///
	ec_exo_polls_times_inc_pp_s ec_pre_polls_times_inc_pp_s ec_end_polls_times_inc_pp_s ///
	using "$OUT\Tables\SHP\ec_all_hhinc_dec_pp_sqrt.rtf", b(3) scalar("N N" "arm1 ABT(1)" "arm2 ABT(1)") se replace
esttab ec_exo_vote_inc_pp_s ec_pre_vote_inc_pp_s ec_end_vote_inc_pp_s ///
	ec_exo_polint_inc_pp_s ec_pre_polint_inc_pp_s ec_end_polint_inc_pp_s ///
	ec_exo_polls_times_inc_pp_s ec_pre_polls_times_inc_pp_s ec_end_polls_times_inc_pp_s ///
	using "$OUT\Tables\SHP\ec_all_hhinc_dec_pp_sqrt.tex", b(3) scalar("N N" "arm1 ABT(1)" "arm2 ABT(1)") se replace label nonumbers booktabs alignment(D{.}{.}{3}) ///
	varlabels ///
				(edu_med "Edu: medium" ///
				edu_high "Edu: high" ///
				age18 "Age" ///
				unemployed "Unemployed" ///
				nlf "Not in Labor Force" ///
				hhinc "Income" ///
				hhinc_dec_pp_sqrt10 "Income (dec.)" ///
				hhinc_sat "Sat. Income" ///
				_cons "Constant" ///
				L.vote "Vote\textsubscript{t-1}" ///
				L2.vote "Vote\textsubscript{t-2}" ///
				L.polint "Pol. Interest\textsubscript{t-1}" ///
				L2.polint "Pol. Interest\textsubscript{t-2}" ///
				L.polls_times "No. polls\textsubscript{t-1}" ///
				L2.polls_times "No. polls\textsubscript{t-2}")





